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We show that a squeezed bath, that acts on the central element of a harmonic chain, can drive the whole 
system to a steady state featuring a series of nested entangled pairs of oscillators that, ideally, covers the whole 
chain regardless of its size. We study how to realize this effect in various physical implementations, including 
optomechanical and superconducting devices, using currently available technologies. In these cases no squeezed 
fields are actually needed, and the squeezed bath is, instead, simulated by quantum reservoir engineering with 
bichromatic drives. 


I. INTRODUCTION 

Quantum technologies and engineering exploit entangle¬ 
ment to boost the efficiency of quantum devices. In spite of the 
many small-scale demonstrations of engineered quantum dy¬ 
namics, the efficient control over large ensembles of quantum 
systems, for realistic realizations of scalable quantum tech¬ 
nology applications, still seems to require significant effort 
from both the experimental and the theoretical side mu. 
Among the many strategies envisaged and pursued for quan¬ 
tum manipulations, reservoir engineering is particularly in¬ 
triguing because it makes use of irreversible/dissipative pro¬ 
cesses to achieve the designed result at the steady state of 
controlled open quantum dynamics These techniques 

have been also applied to the manipulation of arrays, prov¬ 
ing to be in principle applicable to a scalable architecture for 
quantum computation and information {3®. Nevertheless, 
in the schemes proposed so far, many control fields are typi¬ 
cally needed to engineer the dynamics of the many elements, 
and the complexity of the problem increases fast with their 
number, making these strategies of difficult implementation 
with present day technologies. 

On the other hand it has been shown that in certain cases 
it suffices to have access only to a limited part of an en¬ 
semble of interacting elements to control the whole dynam¬ 
ics fl7i M I 22 L l23tl , In this context, here we propose a new 
scheme for the preparation of entangled states in a chain 
of harmonic oscillators with nearest-neighbour coupling, by 
reservoir engineering, with active control over only a sin¬ 
gle element of the chain. Specifically we show that a har¬ 
monic chain can approach an entangled steady state when only 
the central element of the chain is driven by a single-mode 
squeezed reservoir. In particular the oscillators of the chain 
get entangled in pairs. Each pair is formed by the two oscil¬ 
lators at equal distance, on the left and on the right, from the 
central one, as sketched in Fig. [T] Ideally, all the pairs are 
entangled with the same degree of entanglement (which is di¬ 
rectly related to the degree of squeezing of the bath) irrespec¬ 
tive of the size of the chain. This is the result of the interplay 
between localized irreversible dynamics and coherent inter¬ 
actions under specific symmetry conditions. Similar results, 
concerning pairs of entangled systems belonging to two sep¬ 
arated chains have been already discussed in Refs. [17;, 1811 , 


which extended and generalized previous works on the trans¬ 
fer of entanglement to distant quantum objects fl9l - |2T l. How¬ 
ever, we remark that, differently from the scheme discussed in 
Refs, dm where the squeezed bath consists of two en¬ 
tangled modes each of which drives two independent arrays, 
here we consider a single squeezed mode which interacts with 
a single element of the chain. Hence the idea of entanglement 
transfer is absent and, here, no spatially separated entangled 
fields are required. 

We study the performance of this protocol in the case of op¬ 
tomechanical and circuit-QED systems where the harmonic 
chain is realized by, respectively, mechanical and coplanar- 
waveguide resonators. In these implementations no squeezed 
resources are actually needed. In fact, in both cases an ef¬ 
fective squeezed reservoir is realized by coupling the cen¬ 
tral resonator to an additional system driven by a two-tone 
field, which in appropriate parameter regimes can simulate a 
squeezed reservoir in the spirit of reservoir engineering & 
12911 . In the first case the additional system is an optical cavity 
while in the second a superconducting flux qubit. These re¬ 
sults may find application in the context of quantum informa¬ 
tion as efficient strategies to generate entanglement between 
remote, non-directly interacting elements of a quantum pro¬ 
cessor. 

The article is organized as follows. In Sec. [TT] we describe 
the general model and the symmetry of the chains that are re¬ 
quired for the observation of the steady-state entangled pairs. 
The properties of the steady state are then studied in detail in 
Sec. [Till Possible physical implementations are identified in 
Sec.[lv] Finally, Sec.[V]is devoted to outlooks and concluding 
remarks. 


H. GENERAL MODEL: CHAIN OF HARMONIC 
OSCILLATORS WITH A LOCALIZED SQUEEZED 
RESERVOIR 

In this Section we analyze the ideal model that is central 
to this work, and we identify the necessary symmetries the 
model has to fulfill in order to sustain the steady state entan¬ 
gled pairs. We consider a chain of harmonic oscillators with 
nearest-neighbour coupling and with a localized squeezed 
reservoir which affects only the central oscillator. Specifi- 
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FIG. 1: (Color online) An array of linearly coupled harmonic oscil¬ 
lators is locally driven by a squeezed reservoir. The elements in the 
array are labeled by the indices j e [—TV, TV], with j = 0 indicating the 
driven central oscillator. The steady-state entangled pairs are marked 
by dashed arrows. 

cally, we consider a chain of 2N +1 linearly coupled harmonic 
oscillators, with coupling strengths rjj, as depicted in Fig. Q] 
and with creation and annihilation operators b'. and bj, for 
j e {-AT, • • ■ 0, • • -AT), where the index j = 0 corresponds to 
the central oscillator while negative and positive indices indi¬ 
cate, respectively, the oscillators to the left and to the right of 
the central one in such a way that j = -N and j = N indicate, 
respectively, the leftmost and the rightmost oscillator. We as¬ 
sume that the coupling strengths follow a symmetric pattern 
with respect to the central oscillator such that the coupling 
between the oscillators j and j - 1 is equal to that between 
—j and -j + 1. The frequencies of the oscillators, ojj, fol¬ 
low instead an antisymmetric pattern with respect to the cen¬ 
tral oscillator. In detail, pairs of oscillators with indices j and 
—j (for j e {1, • ■ ■ N}) are detuned from the central oscillator 
(the frequency of which is u>q) by opposite quantities, that is 
tOj - too — a>o - a)-j = A j. The effects of deviations from 
this symmetric situation are discussed in the next section. The 
corresponding Hamiltonian, in the interaction picture with re¬ 
spect to H 0 = h u> o Ey=-w b'j bj , is given by 

N 

H cham = hj^ A j (b) bj - b\. b-j) (1) 

7=1 

N 

+n 2 r\j(b h 1 b) + b) , bj + b. j+1 b'' . + // /+1 h /) . 

7=i 

As shown in detail in the next section, by virtue of these 
symmetry conditions the chain is invariant under the effect of 
specific Bogoliubov transformations. The corresponding Bo- 
goliubov modes are combinations of pairs of oscillators, and 
are the modes that are actually stabilized by the effect of the 
squeezed bath, that drives, with rate F, only the central oscilla¬ 
tor. Its effect on the system density matrix p can be described 
in terms of a Lindblad operator, of the form [30] 

4 r 4£p = r {(n + 1) £> [*o, bl] + n <D [fcj, b 0 \ 

-mD[bo,b 0 ]-mD[bl,bl]}p (2) 

with 

T)\x,y\p = 2 xpy —y xp -py x . (3) 


This model describes, for example, the effect of the output 
field of a degenerate parametric oscillator (similar to the sys¬ 
tem studied in Ref. [ l_8j]). The statistical properties of the 
squeezed reservoir are determined by the parameters n, which 
accounts for the number of excitations, and in, which ac¬ 
counts for the correlations. In general |w| is bounded by the 
value |w| < x//T(/? + 1), with the equality corresponding to a 
pure state of the reservoir. When \m\ < \jh{n + 1), instead, 
the reservoir is in a squeezed thermal state. Specifically, the 
quadratures of the reservoir are squeezed below the vacuum 
level only when |m| > n, and the variance of the maximum 
squeezed quadrature (relative to the vacuum noise level) is 

S = 2 n + 1 - 2 |m|. (4) 

Finally in our description we also take into account thermal 
dissipation at rates jj in a bath with nj t thermal excitations, 
which can be described by a Lindblad operator of the form 

N 

LdissP = Yj Vi \{ nr J + l ) !D [ b i- b )] + ,lT J n i b )’ b j]} p ■ 

j=-N 

(5) 

In summary we study a system described by the following 
master equation 

P = \H ch ain,p\ + + £di«P ■ ( 6 ) 

We will show that, at the steady state, oscillators at equal dis¬ 
tance, on the left and on the right of the central one, are entan¬ 
gled in pairs as depicted in Fig. [I] 

We remark that the squeezed bath, described by Eq. Q, 
is equal to a standard dissipative bath for a Bogoliubov 
(squeezed) mode with annihilation operator c r foo-e 1( °s,-fcy, that 
is 

£%tfzP = r {(«r +1)2) \c r bo - fPSrbl ^ S r b {] \ 

+ii r D [c r b f 0 - e" 1(C i,. b 0 , c r b 0 - e v .v4]}p - (7) 
with e v = in/ \rh\, s r /c r = (n - h r ) / |m|, c 2 r - s 2 r = 1 and 



Correspondingly, as a result of the symmetric pattern of cou¬ 
plings jjj that we have chosen, the central Bogoliubov mode 
is coupled to the two modes c r b ± \ + e I<J s r 7>^ l and in turn these 
are coupled to similar modes realizing a chain of Bogoliubov 
modes whose annihilation operators are c r b ±j - (-l)- , e 1 ^5,i)i.. 
Furthermore we note that the Hamiltonian term proportional 
to the detunings A,, in Eq. (Q1 does not induce additional cou¬ 
plings between these Bogoliubov modes. This is a conse¬ 
quence of the chosen antisymmetric structure of the detunings, 
and of the fact that, in general the number difference operator 
b\bj - b^bk is invariant under the effect of a two-mode squeez¬ 
ing transformation. In conclusion, as discussed in more detail 
in the next Section, due to the specific symmetries that we 
have selected, the squeezed bath acts as a standard dissipative 
reservoir over a chain of Bogoliubov (squeezed) modes. 
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III. STEADY STATE ENTANGLEMENT 


A simple analytical solution for the steady state is found 
when -Cdiss = 0 (i.e., jj — 0 V;). In this case it is useful to 
study the system in a new representation defined by the unitary 
operator 


U = e (-lV>(en] bl^bj b-j) t 

which realizes the Bogoliubov transformation 

U bo I/ 1 = cosh (r) bo + e"* 5 sinh (r) 

UbjU 1 = cosh (r) bj + (-1 sinh (r) b^_ ., 


(9) 


( 10 ) 


with tanh(r) = (n - h r ) /\m\ and h r defined in Eq. We note 
that h r = 0 when \m\ = y/h(h + 1), that is, when the reservoir 
is in a pure squeezed state. 

In the new representation the system Hamiltonian remains 
unchanged as a result of the symmetry that we have chosen. 
The master equation for the transformed density matrix p = 
U'p U is, therefore, given by 


P — , [H chain , p\ 

n 

+r {(hr + 1)D [bo, b\\ + h/D [bl bo\}p (11) 

that has the same structure of Eq. 0 but with m = 0, and h r in 
place of h. This master equation describes a system interact¬ 
ing with a thermal environment with h r thermal excitations. 
In this representation, at large times, the system thermalizes 
and all the oscillators approach the same thermal state with 
h r excitations p^ y . = £„ j±r- (^-) |«)<w| (equal for all ;'), so 
that the global steady state is given by 

Pst = dj=- N P% > ( 12 ) 

which is equal to the vacuum when h r = 0, that is when the 
state of the driving reservoir is pure. We have to remark, 
however, that although this thermal state is always a steady 
solution, in certain situations it is not unique. If for exam¬ 
ple the matrix of coefficients, AI, of the system of equations, 
P(t) = A\P(t), for the evolution of the average field operators 
P = ((b - N ), • • • , (bo) , • • • , (bn)) 7 , corresponding to Eq. (fTIT) . 
has eigenvalues with zero real part, then the corresponding 
eigenmodes do not dissipate, and the steady-state solution is 
not unique. This happens, for example, when there exists an 
eigenvector, v*, of the matrix 


-A n r/ N 

0 



0 1 


0 

m 

0 

m 

0 

(13) 

0 


0 

r/N 

m 

A n , 



that is orthogonal to the vector 

v 0 = (0, ■ ■ ■ 1, ■ • -0) T (14) 


that corresponds to the central oscillator, i.e., vj • vo = 0. In 
this case, in fact, the corresponding eigenmode remains unaf¬ 
fected by the thermal reservoir. A specific realization of this 
situation is observed, for example, when all the oscillators are 
resonant A j = 0. If instead all the eigenvectors, v*, of Al|^ ^ 

are non-orthogonal to vo, i.e., vj ■ Vo 4 0, then the steady state 
identified above is unique. We have checked that this condi¬ 
tion is fulfilled when, for example, the couplings are uniform 
and the detunings vary linearly according to the relations 

Vj = '?> for; e {1 ••■A) , 

A j = A + jS, for j e {l ■■■ N] . (15) 

We anticipate that the specific results that we discuss hereafter 
are obtained for chains which follow this structure. 

Before studying the system dynamics under more general 
conditions we have to analyze the entanglement properties of 
the steady state that we have just identified, and that is valid 
when -Cdiss = 0. Specifically, in the original representation, 
the thermal steady state in Eq. (fl2l >. reduces to p st = U p st U\ 
with U defined in Eq. (0. that corresponds to a factorized 
state of the form p st = ®j =0 Pst,j with the central oscillator in a 
squeezed thermal state 

Pst,o = Up%U f , (16) 

and all the pairs of oscillators with opposite indices, j and 
in two-mode squeezed thermal states 

Pstj = Up% (17) 

that is the same for all the j regardless of the size of the chain. 
We note that these density matrices describe pure two-mode 
squeezed states when h r = 0. The corresponding logarith¬ 
mic negativity I3ll l32j between each pair of oscillators can be 
expressed in terms of the corresponding correlation functions 
(see App.[0 and is given by 

En[ j, ~j] = max {0, - log 2 S} , (18) 

where S is defined in Eq. 0, which shows that the strength 
of the pairwise entanglement increases with the degree of 
squeezing of the driving reservoir. We remark that this value 
is independent from the size of the chain and from the param¬ 
eters rjj and A j. 

The results that we have discussed so far have been found 
for -Cdiss = 0. When this condition is not satisfied it is, in gen¬ 
eral, not possible to identify a simple analytical solution. In 
this case we resort to the numerical evaluation of the steady 
state and we characterize its entanglement properties in terms 
of the logarithmic negativity 113 lLl32l l for different pairs of os¬ 
cillators £, v [;j, ji]- In general finite values of yj tend to re¬ 
duce the entanglement of the pairs, which, however, remain 
visible as long as the thermal noise is not too strong, as de¬ 
picted in Fig. [2] (a) and (b). In any case, only oscillators 
with opposite indices can be entangled, i.e. E^[ji,j 2 ] = 0 
for ;j 4 -ji- We further note that, when yj 4 0, the pair¬ 
wise entanglement depends non-trivially on the structure of 
the chain. For example, when the chain structure is described 
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FIG. 2: (Color online) (a): Steady-state logarithmic negativity, E.vly'i , ji\, for all the pairs of oscillators in a chain of 9 oscillators (N = 4), 
with couplings and detunings defined in Eq. at the steady state of the dynamics described by Eq. with n = 2, m = Vn(n + 1) 
(corresponding to reservoir squeezing of ~ 10dB), 77 = T, A = -0.4 T, 6 = 0.1 T, and uniform dissipation with y ; - = 0.01 T and n T j = 0 V j. 
(b), (c), and (d): Steady-state logarithmic negativity for the pairs with indices {j, -j}, as a function of, respectively, the mechanical dissipation 
y = jj, V j, the oscillators couplings 77 , and of the detuning A. The other parameters are as in plot (a). The vertical thin lines in (b), (c) and 
(d) indicate the parameters corresponding to plot (a). The horizontal thin lines indicate the level of entanglement at yj = 0, Vy, as defined in 
Eq. J 1 8b . In all cases the other pairs are not entangled, i.e. E N [yj, y' 2 ] = 0 for j\ —ji. 



tj (units of T) 


A (units of O 


k=9 

k=8 

k=7 

k=6 

k=5 

k=4 

k=3 

k=2 

k=l 






N 

A 

— 

1 

-r 

— 

2 

-0.57r 


3 

-0.48r 


4 

1 

0 

Ji. 

— 

5 

- 0 . 42r 

— 

6 

-0.43r 

— 

7 

- 0 .46r 

— 

8 

-0.495r 

— 

9 

-d.54r 


FIG. 3: (Color online) Projection • Vo of the eigenmodes of the 
chain, described by the eigenvectors v* of the matrix Al| r , defined 
in Eq i fTTl . over the central oscillator, described by the vector Vo, 
defined in Eq. 03- as a function of (a) the oscillators coupling 77 , 
and (b) the detuning A. Increasing values of the index k correspond to 
increasing values of the frequency of the normal modes. The vertical 
lines indicate the optimal parameters where all the curves converge 
to a close value, and correspond to the parameters of Fig. [2] (a) and 
to the vertical lines in Figs. [2] (c) and (d). The other parameters are 
as in Fig. 0(a). 


FIG. 4: (Color online) Steady-state logarithmic negativity, 
for the pairs of oscillators with indices j and — j, for chains of in¬ 
creasing size (from 3 to 19 oscillators, i.e. (V e {1 ■ ■ • 9)), as detailed 
in the table. For each size the detuning A has been optimized by 
matching the configuration for which the projections v* • Vo are as 
close as possible according to the procedure outlined in Fig. [3] For 
each j the lowest mark corresponds to the entanglement realized with 
the largest size {N = 9). Increasing values of E n correspond to de¬ 
creasing sizes. The dashed lines connect the results obtained with the 
same value of N. The other parameters are as in Fig. 0(a). 


by Eq. (IT5] >. the entanglement of the pairs is maximized for fi¬ 
nite values of 77 , while it oscillates as a function of A (see Fig. [2] 
(c) and (d)). In detail, the optimal value of the entanglement is 
observed when the projection \ T k ■ Vo of the eigenmodes of the 
chain, vy-, over the central oscillator, Vo, are of similar size as 
described by Fig. [3 In this case, in fact, all the normal modes 
are equally coupled to the reservoir, and they are all efficiently 
driven to their steady states. The behaviour of the steady state 
as a function of the size of the chain is shown in Fig. [4] As 
already stated when yj = 0 all the pairs share the same degree 
of entanglement regardless of the size of the chain. In con¬ 
trast, as shown in Fig. [4] when -C,/ !SS + 0, the entanglement of 
each pair decreases mildly with the number of the oscillators 


in the chain. We also observe that, similar to the findings of 
Ref. HI, the entanglement decreases with the distance of the 
pair from the central oscillator, and exhibits a weak revival for 
a few pairs at the ends of the chain. 

We finally discuss the effects of random deviations from the 
symmetric configuration identified in Sec. QI] The outcome 
of this analysis is reported in Fig. [3 where we show the re¬ 
sults for the steady state entanglement evaluated considering 
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FIG. 5: (Color online) Steady-state logarithmic negativity, 

E n [j, - for the pairs of oscillators with indices j and —j. Each 
plot reports the results for 200 different realizations (red-solid- 
thin marks) of the dimensionaless random variables £ n j and 
uniformly distributed in the ranges £n,j 6 and j s 

[ - <5<Ta,<5£a], which describe small fluctuations of the oscillators cou¬ 
plings and detunings about the ideal symmetric situation identified 
in Sec. mi as discussed in the main text. Each plot is evaluated with 
different values of the range of variability 8 £, and 8fa of the ran¬ 
dom variables, as specified in the figure. The black-dashed-thick 
marks indicate the entanglement in the symmetric configuration, i.e. 
£, h j = £ A j = 0, V j. The other parameters are as in Fig.[2](a). 


a Hamiltonian of the form 

N 

H Chain = h ^ (A + jS) [(l + j) b\ bj - (l + £ A ,-j) b\j b ; | 
j= 1 

N 

+* 2>[( 1 * £*-j) & 1 b ) +h U b j) 

j= 1 

+{'+^-j){b- J+t bi j+ b^b-j)\ , d9) 

where (, hJ and £a j are random variables uniformly distributed 
in the ranges <5^] and 1 5£ A ] respectively. Specif¬ 

ically, Eq. ( fl9l > accounts for deviations from the configuration 
described by Eq. (Q} with the parameters defined in Eq. (ITfil i. 
Fig- 0 demonstrates that the protocol is significantly resilient 
to small asymmetries. We also note that the protocol is much 
less sensitive to deviations in the detunings than in the oscil¬ 
lators couplings. 



FIG. 6 : (Color online) Steady-state logarithmic negativity E N [j,-j] 
for the pairs of resonators with opposite indices, in a chain of 9 (N = 
4) mechanical resonators, as a function of (a) |G_|, with k = 0.01w 0 , 
and (b) k, with |G_| = 0.85 |G+| (£_ = 136wo). The solid and dashed 
lines, that are almost indistinguishable because very close to each 
other, are evaluated from the steady state of, respectively, Eq. d 
and Eq. (( 6 } with the parameters defined in Eq. <221 . The vertical line 
in (a) indicates the parameters used in (b) and similarly the one in (b) 
indicates the value of k used in (a). The other parameters are ato = 
2jtx 1GHz, g = 5x 10~ 5 «o, 8 = 0.2x 10~ 3 tno, A cav = 1550nm, |G+| = 
8 x 10 -3 Wo (£ + = 160 aj 0 ), temperature T = 50mK and y, = 10~ 5 to 0 
y j. Moreover the chain follows the structure defined by Eq. G3 
with r] = 2 x 10~ 3 a>o and A = 0.8 x 10~ 3 Wo that are consistent with 
the optimal values identified in Fig. [3] All the other pairs are not 
entangled, i.e. E N [j u j 2 l = 0 for j, * -j 2 . 


parametric oscillator below threshold, which would consti¬ 
tute the squeezed reservoir similar to the proposal discussed 
in Ref. [18]. A promising alternative approach, that we ex¬ 
plore here, is instead that of realizing a squeezed bath by reser¬ 
voir engineering using bichromatic driving fields. The idea of 
using bichromatic fields to engineer quantum dissipative dy¬ 
namics has been introduced in the context of trapped ions in 
Ref. @]. It has been then extended to a number of different 
physical setups including mechanical resonators j25l - [28il and 
circuit-QED 112911 . All these systems are in principle suitable 
candidates for the engineering of the dynamics that is the sub¬ 
ject of the present proposal. Hereafter we briefly review and 
extend these techniques based on two-frequency drives to the 
simulations of our protocols with two specific setups. First, 
we consider an array of mechanical resonators and then one 
of superconducting microwave resonators. In both cases we 
identify the parameters for which the dynamics discussed in 
Sec. mi can be efficiently simulated. We characterize the sys¬ 
tem dynamics in terms of the logarithmic negativity for pairs 
of oscillators evaluated at the steady state of the system dy¬ 
namics showing that the results identified in the previous sec¬ 
tion can be reproduced under specific parameter regimes. 


IV. PHYSICAL IMPLEMENTATION: RESERVOIR 
ENGINEERING WITH BICHROMATIC DRIVING FIELDS 

Here we analyze how to implement the dynamics discussed 
in the previous section with actual physical systems. We 
first note that, in principle, one could realize the model de¬ 
scribed in Sec. [TT]by, for example, driving arrays of optical or 
microwave resonators with the output field of a degenerate 


A. Optomechanics: Simulation with an array of mechanical 
resonators 

In this section we describe how to engineer the dynamics 
introduced in Sec. uni with a chain of mechanical resonators 
similar to those discussed in |[33M35ll , which are based on 
optomechanical crystal implementations i36U4l!l . where me¬ 
chanical and optical modes are defined by localized defect in 
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and the mechanical resonators about their corresponding av¬ 
erage fields, and the linearized coupling strengths G+ and G 
are proportional to the intensity of the driving fields on the red 
and blue sideband respectively (their specific dependence on 
the other system parameters is discussed in Ann. IB 1 b . 

As discussed in [27}], when |G+| > |G_| (which guarantees 
the stability of this model, keeping however in mind that in 
certain cases, when |G+| ~ |G_|, non resonant processes that 
have been neglected in this effective model (see App.|B} can 
modify the system stability l42ll ). this model describes the 
cooling of the mechanical Bogoliubov mode of the central os¬ 
cillator given by Bq = (G + ba + G-b^j / V|G >| 2 - |G_| 2 . If the 
cooling dynamics is sufficiently fast to overcome the thermal 
noise then such a mode can be cooled to its ground state, that 
corresponds to a squeezed state for the original mode. In turn 
the other oscillators of the chain are driven to an entangled 
steady state following a dynamics similar to the one discussed 
in Sec. [Ill] In particular, if the decay rate of the cavity field, 
k , is sufficiently large, the degrees of freedom of the optical 
cavity can be adiabatically eliminated and the system can be 
described by a model equal to Eq. © with 


FIG. 7: (Color online) Steady-state logarithmic negativity EnU’-J] 
for pairs of resonators with opposite indices, in a chain of 9 mechan¬ 
ical resonators as a function of the dissipation rate of the resonators 
(assumed uniform through the chain, yj = y V j), for three different 
values of the temperature T as specified by the arrows. For each set 
of points the linearized coupling G has been chosen in order to op¬ 
timize the entanglement of the pairs, and the corresponding values 
of £ = |G_| / |G + | are indicated in the plot. The solid and dashed 
lines connect the results found from the steady state of, respectively, 
Eq. ( 120b and Eq. ©, where in the latter the bath parameters are de¬ 
fined in Eq. ( 122b . The other parameters are as in Fig. [6] All the other 
pairs are not entangled, i.e. E,vlii > 721 = 0 for j\ t- -ji- 


planar artificial quasi-periodic nanostructured crystals. The 
system hence consists of an array of localized mechanical 
modes in a phononic-crystal. Phononic excitations can hop 
between nearby mechanical modes [35] realizing the chain. 
The squeezed reservoir is, instead, realized by having the cen¬ 
tral mechanical mode coupled to an optical mode iPioi. l4lj| 
that, in turn, is driven by a two-tone field realizing a reservoir¬ 
engineering scheme similar to the one investigated in | f26li27ll . 
The corresponding linearized optomechanical dynamics can 
be described by a master equation for the system density ma¬ 
trix p of the form (see App.[B]for a detailed derivation) 

p = -^\H,p\+(£ a +£di SS )p, (20) 

Ti 

where, the cavity dissipation at rate k is described by £„ = 
kD [a, a^J, with D defined in Eq. ©, -Cdiss is introduced in 
Eq. ©, and the Hamiltonian is given by 

H = H chain + n [a t (G + b 0 + G-bl ) + a (G* + ^ + G*_b 0 )] .(21) 

Here H c i win is defined as in Eq. (jTJ, the operators a and /?, ac¬ 
count for the fluctuations of, respectively, the optical cavity 


r 


n 


m 


IG+I 2 ~ IG-I 2 

K 

\G-\ 2 

\G+\ 2 - |G _| 2 

g; g_ 

K ? + | 2 - |G _[ 2 ' 


( 22 ) 


Here we consider systems similar to the ones investigated 
in iMS where the mechanical modes are at gigahertz fre¬ 
quencies while the optical ones are in the near infrared. In 
detail we choose ljq = 2nx 1GHz and the wavelength for the 
cavity field A/ = 1550nm. Figs.[6]and[7]show that the model 
of Sec. [IT] describes with high accuracy the dynamics of this 
system in the corresponding parameter regime. We observe, in 
fact, that dashed and solid lines, corresponding to the results 
obtained from the steady state of, respectively, Eqs. © and 
m, are always very close. The pairwise entangled steady- 
state is achieved whenever the effective squeezed reservoir 
is efficiently realized. The parameter regimes for which the 
bichromatic drives squeeze efficiently a mechanical oscillator 
have been discussed in detail in Ref. j27ll . In particular, the 
squeezing is optimized at specific values of |G_|/|G+[, and 
of the decay rate of the cavity field. Correspondingly we ob¬ 
serve strong pairwise entanglement as described by Figs[ 6 ](a) 
and (b). In turn, when the optimal conditions for G± and k are 
satisfied, the entanglement is achieved if the chain follows the 
symmetries identified in Sec. [II] and the couplings and detun¬ 
ings are optimized as discussed in the previous section. In 
this case, all the pairs of resonators with indices j and —j 
exhibit roughly the same value of the logarithmic negativity 
£,v [j, —j], and in general larger entanglement is achieved for 
larger quality factor u>o/yj of the resonators and for smaller 
temperatures as shown in Fig. [7] 
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FIG. 8: (Color online) Steady-state logarithmic negativity 
for pairs of resonators with opposite indices, in a chain of 9 copla- 
nar waveguide resonators as a function of the dissipation rate of the 
resonators (assumed uniform through the chain, = y V /'), for three 
different values of the temperature T as specified by the arrows. The 
coupling strength G_ has been chosen in order to optimize the corre¬ 
sponding entanglement for each set of points, and the corresponding 
values of £ = |G_| / |G + | are indicated in the plot. The results are 
evaluated from the steady state of Eq. @ with the bath parameters 
defined in Eq. 122b . and with the effective coupling G ± defined in 
Eq. 126b . The other parameters are coo = 2n x 5GHz, g = 0.06wo, 
r/ = 2 x lCT 3 Wo. A = 0.8 x 10~ 3 w 0 , 6 = 0.2 x 10 -3 «o, e = 2a>o, 
|G + | = 6 x 10~ 3 Wo (£+ = 0.25tuo) k = 0.02a>o. All the other pairs are 
not entangled, i.e. E N [j u j 2 ] = 0 for j, * -j 2 . 


B. Circuit QED: Simulation with an array of superconducting 
coplanar waveguide resonators 


A different system that can be used to engineer the dynam¬ 
ics described in Secs. QI] and |III]consists of a chain of capac- 
itively coupled microwave coplanar waveguide (CPW) res¬ 
onators B431144 1. This system is attractive because CPW res¬ 
onators can be efficiently coupled to superconducting qubits, 
hence realizing a dynamics similar to the one discussed in 
Ref. |d, where also the qubits get entangled in the steady 
state. 

In this case, the squeezed reservoir can be engineered 
by coupling only the central resonator to a superconducting 
qubit. Specifically one can implement the scheme described 
in Ref. l29ll . where a superconducting flux qubit |45. 46| with 
tunable gap j47l ]48il . is coupled to a CPW resonator at fre¬ 
quency a>o, and the gap is modulated by a bichromatic field. 
We remark that in principle similar dynamics can be realized 
with different kinds of superconducting qubits. At the degen¬ 
eracy point the flux qubit is described (following the notation 


in Ref. f29il ) as an effective spin 1/2 in a modulated magnetic 
field 


H q = n 


— - fi+ cos(e + u>o) - & cos(e - wq) 


cg 


(23) 


while the interaction with the resonator, the creation and an¬ 
nihilation operators of which are bo and b' Q , is given by 

Hi - Tig (o"+ + o'-) (bo + bl) , (24) 

where cr ± are the lowering and rising operators for the effec¬ 
tive spin. Here we assume that this resonator is the central 
element of a chain of CPW resonators described by a Hamilto¬ 
nian equal to Eq. (IB 111 . In the interaction picture with respect 
to H q + h (Jo Yjj bjbj one finds that, similar to the derivation 
in Ref. |2Sj], the total system Hamiltonian at lowest order in 
£ ± /(e ± wo), and retaining only the dominant resonant terms, 
under the condition k, gS ± /(e ± wo) <sc tJo, reduces to 

H s = H c h a i„ + Ti(r + [ G+b 0 + G-Z?, 1 ,] + h.c. (25) 

where 


G ± = , (26) 

e±a>o 

and H c hai„ is defined as in Eq. {Q>. Also in this case, if the 
qubit decays fast enough, then the Bogoliubov mode Bo = 
(G+bo + G-b J) / V|G+| 2 - |G_| 2 can be efficiently cooled. Fast 
relaxation, that should overcome the natural dephasing of the 
device, can be engineered by constructing the sample with a 
nearby transmission line which serves as an additional dissi¬ 
pation channel. If the resulting relaxation rate, that we denote 
as k, is sufficiently large then the qubit remains mostly in its 
ground state. In this case the effective spin can be approxi¬ 
mated as a harmonic oscillator. Thereby the system descrip¬ 
tion reduces to a model equal to Eq. ( l20l) . where the bosonic 
operators a and a 1 have replaced cr_ and cr+, respectively. We 
have verified that using this model, with the parameters that 
we will discuss below, the population of the bosonized qubit 
is always very small (< 0.02), hence providing a justification 
of the results presented hereafter. Specifically, under these 
conditions, the degrees of freedom of the qubit can be adiabat- 
ically eliminated, and also in this case the system is described 
by a model equal to Eq. ©, with the bath parameters defined 
as in Eq. (l22li . 

In detail, we consider system parameters consistent with 
state-of-the-art circuit-QED experiments |44jj. We focus on 
resonators with frequencies of u>o — 2n x 5GHz, and qual¬ 
ity factors Q ~ 10 5 . The frequency of the central qubit 
is 2n x 10GHz and its coupling to the resonator is g ~ 
2n x 300MHz. Moreover we employ a large decay rate of 
the qubit k ~ 2n x 100MHz, much larger than its natural de¬ 
phasing rate (the largest dephasing time in flux qubits is of 
~ 10/zs m, corresponding to rates of tens of kilohertz) which 
here is therefore neglected. Assuming a chain of resonators 
described by coupling strengths and detunings consistent with 
those identified in Sec. [Till the steady-state results are very 
similar to those found for the mechanical resonators as de¬ 
picted in Fig. [8] 
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V. DISCUSSION AND CONCLUSIONS 

We have shown that a harmonic chain with nearest- 
neighbour coupling can be driven to a non-trivial steady state 
exhibiting a series of nested entangled pairs covering the 
whole chain when only the central oscillator is coupled to a 
single mode squeezed bath. Ideally, when the efFects of ad¬ 
ditional noise are negligible (jj ~ 0 V j), the entanglement of 
all the pairs is the same regardless of the size of the chain, 
and its strength is directly related to the degree of squeez¬ 
ing of the reservoir (see Eq. (IT8]>). We have shown that the 
efficiency for the production of the pairwise entanglement is 
strongly dependent on the structure of the chain. In particular 
we have identified the symmetry conditions which are nec¬ 
essary for the observation of this phenomenon, and we have 
described, in a simple specific realization, how to optimize the 
entanglement based on the study of the normal modes of the 
chain. Finally, we have studied how to simulate this effect, 
using bichromatic drives, and with actual physical systems, 
hence demonstrating that optomechanical as well as circuit- 
QED systems offer the potentiality to harness the discussed 
dynamics in parameter regimes that are within the reach of 
near-future experiments. 

A. Alternative implementations 

Let us now comment on the possibility to simulate the dy¬ 
namics described in this article with other physical systems. 
We remark that any system that can be mapped, in some limit 
to a harmonic chain is a suitable candidate. Furthermore, we 
observe that in principle the squeezed reservoir can be realized 
using the output held of a degenerate parametric oscillator op¬ 
erating below threshold (similar dynamics have been studied 
for example in Ref. ifTiilP . It consists of a continuous-wave 
squeezed held mode that, hence, exhibits entanglement be¬ 
tween the spectral components at opposite sideband frequen¬ 
cies [50]. In this conhguration, the entangled pairs production 
that we have analyzed can be seen as a means to extract and 
spatially separate the entanglement contained in a single mode 
squeezed continuous held. 

Other systems which can be manipulated with bichromatic 
drives, and that can be exploited in our protocol are listed be¬ 
low. 

a. Array of optical cavities with optomechanical manip¬ 
ulation: First we note that in principle one can envisage 
an implementation, complementary to the one studied in 
Sec. IIV Al which makes use of an array of optical cavities 
where only the central one interacts with a mechanical res¬ 
onator. The central cavity is driven by a two-tone held to en¬ 
gineer the squeezed bath similar to the results discussed in 
Ref. |2 81. The system is thus described, with minor modifi¬ 
cations, by the model of Sec. IIV Al (see also App. IB 21) . In 
this case, the chain of optical cavities can also be realized 
with arrays of microcavities which are coupled via evanescent 
waves l5lll or by optical-hbers (52]. 

b. Arrays of mechanical resonators: Apart from op¬ 
tomechanical crystals, other realization of mechanical arrays 


based on nano- and micro-mechanical cantilevers, beams, and 
membranes, and that can be manipulated with either opti¬ 
cal, electrical, or magnetic forces, could be used for our pur¬ 
pose ll53l - l58ll . In particular electromechanical arrays are very 
promising systems that are actively investigated in many ex¬ 
periments. They can be efficiently coupled to superconduct¬ 
ing microwave cavities [!59|, 60] which can, thus, be used 
to realize a model formally equal to the one described in 
Sec. IIV Al Moreover squeezed reservoirs for an electrome¬ 
chanical resonator have been recently realized using two- 
frequency drives min. 

c. Array of superconducting microwave resonators with 
Electro-optical manipulation: Another intriguing possible 
implementation relies on the realization of the squeezed reser¬ 
voir, for a chain of CPW resonators, by coupling the cen¬ 
tral oscillator to an optical cavity via an electro-optic mod¬ 
ulator |j 6 l|. It consists of a second-order nonlinear-optical 
medium which is simultaneously placed inside an optical cav¬ 
ity and between the plate of a capacitor which is connected 
to the CPW resonator. It induces a voltage-dependent phase 
shift to the cavity optical held, whose corresponding interac¬ 
tion Hamiltonian is, thus, equal to the standard optomechan¬ 
ical one. Consequently, if the cavity is driven by a two-tone 
held, this system is described by a model equal to the one 
discussed in Sec. IIV Al and in App.[B] The experimental fea¬ 
sibility of similar devices is still to be demonstrated, and the 
actual implementation of these ideas would, most probably, 
require careful studies of improved system design in order to 
achieve sufficiently large coupling strengths. 

d. Trapped ions: Trapped ions are highly controllable 
systems that have been proved to be an efficient platform 
for quantum simulations and engineering |{ 62 jl . and for which 
the squeezed reservoir, engineered with two-frequency drives, 
has been already experimentally realized [10]. In this case, 
however, the challenging aspect is probably the design of the 
chain. In fact, the long range nature of the Coulomb inter¬ 
action between the ions makes the realization of chains with 
only nearest neighbour coupling problematic. In any case, it 
is worth observing that, this aspect does not necessarily pre¬ 
vent the possibility of applying our protocol to trapped ions. 
Indeed, an interesting question, that deserves further investi¬ 
gation, is how the dynamics that we have studied is affected 
when long range interactions are taken into account. 


B. Outlook 

In conclusion we observe that, while we have shown how 
to optimize the entanglement of the pairs under the specific 
simple conditions defined by Eq. 0 , a still open and com¬ 
pelling question is the identification of the optimal structure 
of the chain (coupling strengths and frequencies of the oscil¬ 
lators) which could result in the maximum steady-state entan¬ 
glement depending on the specific statistical properties of the 
driving reservoir. Moreover, we note that it would be inter¬ 
esting, in relation to the analysis of Sec. IIV Bl to extend the 
presented scheme to the situation in which each resonator in¬ 
teracts with a superconducting qubit, such that also pairs of 
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qubits could get entangled at the steady state, analogously to 
the findings discussed in lfl7ll . Similar considerations in the 
context of trapped ions could, instead, lead to the creation of 
entanglement between the internal states of the ions. Finally, 
another very challenging, but fascinating, future research di¬ 
rection would be the investigation of possible extensions of 
this protocol for the generation of more complex quantum 
states, in higher dimensional arrays and with arbitrary-range 
interactions which could be useful, for example, in the context 
of one-way quantum computation. 
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tonian of the form 


Hchain — bl b^b{) + tl EE + fAy) 

f=±l 3 = 1 




+h X X 7,1 1 b l(j~ i) + h f (j ~ 11 11 b f J + b £i l (B ' 

f=±l j= 1 


where positive and negative indices indicate resonators that 
are placed, respectively, to the left and to the right about the 
central one whose index is j = 0. The central oscillator is 
coupled to an optical cavity at frequency u> c with annihilation 
and creation operators a and a\ via an optomechanical inter¬ 
action, with strength g , described by 

H om — ftga^a (b^ + b 0 J ■ (B2) 

The cavity is driven by two laser fields at the mechanical 
sideband frequencies o>l + too, where the central frequency 
is shifted from the cavity resonance by a small quantity 


Appendix A: Steady state correlations when £ diss = 0 


6 C = CO c - (1>L , 


Here we study the correlations between the operators of the 
oscillators at the steady state of the dynamics described by 
Eq. © when Jidiss = 0. Since the system is Gaussian, the 
following expressions fully characterize the steady state. 

Specifically, the correlation matrix for the central oscillator, 
whose elements are written in terms of the vector of operators 
b 0 = (bo,b ] ^ as {C 0 } fr = Tr[{b 0 }^ (b 0 A p st , 0 ], where p. s ,, 0 
is defined in Eq. (IT6l > and E, V are vector indices not to be 
confused with the indices of the oscillators, is 


Co 


m n + 1 
n in* 


(Al) 


This corresponds to a squeezed state if \m\ > n, that is when 
the reservoir is actually squeezed. The correlation matrix 
for two modes j and — j, the elements of which are written 
in terms of the vector of operators by = (bj, b-j, b 7 , b' j as 


(<4, = T '- 

(N, M, 

Pvjj, 

with p stJ 

defined in Eq. (IT71). 

is, instead given by 






° (- 

-1 )hn 

n + 1 

0 ' 


Cj = 

(— 1 )hn 
n 

0 

0 

0 

0 

n + 1 

(-1 ym* 

(A2) 


0 

n 

(-1 ym* 

0 J 



that corresponds to an entangled state when \m\ > n, and de¬ 
scribes a pure state when \m\ — V n(n + 1). All the other cor¬ 
relation functions between different oscillators are zero. 


Appendix B: Optomechanical chain driven by a bichromatic 
field 


We consider 2 N + 1 resonators, in a configuration similar 
to the one discussed in section HI1 and described by a Hamil- 


the specific value of which, as specified below, is chosen to be 
opposite to the shift of the cavity resonance due to the optome¬ 
chanical interaction. The corresponding Hamiltonian term de¬ 
scribing the cavity driving, in a reference frame rotating at the 
frequency col, is given by 

H d = ft a 1 e~ i<Jo ' + £ + e i<J °'] + h.c ., (B3) 

where h.c. indicates the Hermitian conjugates. We further 
take into account the dissipation of the cavity field at rate k, 
which is described by the standard Lindblad operator £„ = 
k'D |t/,a + |, where D is defined in Eq. <[3j, and the thermal 
noise at rate jj, due to thermal reservoirs with njj thermal 
excitations, that is described by Eq. ©. 

The total system dynamics is therefore described by the 
master equation 

P — —— + Hchain H om + ! 1 d , p J + {Jba + Jldiss) P ■ 

(B4) 


1. Linearized dynamics, approximation in powers of g and 
elimination of the fast rotating terms 

The optomechanical non-linear dynamics can be linearized 
around the steady state average amplitudes of the cavity field, 
a(t ), and of the mechanical resonators, fij(t), which fulfill the 
relations 


a(t) = - (i 6 C + k) a(t ) - 2 i g a(f)Re (A)(0) 

-i fs_e' iw "' +£ + e'"“'] 

J3(t) = - (HV + i cu 0 - i g |ar(0l v o > (B5) 

where we have introduced the vectors 

m - os-,v(o, • • •#>(*), • ■ -MO) 7 ', 

vo = (0, - - -1, - ■ ■ 0) r , (B6) 
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and the matrix 


w- N i ij N 0 


0 \ 


TV = 


iriN 

0 


ir/i w 0 i /71 0 




0 


'' • i r/N 

0 \rj N w N , 


(B7) 


with w±j = +iA y - + y ± j for j e {0,1, • • -N}. In the shifted rep¬ 
resentation described by the density matrix p' - U(fp p U(t), 
where U(t) is the displacement operator U(t) = D a ( t ) 

D m , with D a(t) = e “« ot -“®* 0 and D m = ePiW b }-W b *, 

that realizes the transformations 


U(t)^aU(t) — a + a(t) 

U(t) f bj U(t) = bj+fijit), (B 8 ) 

the optomechanical interaction Hamiltonian can be approxi¬ 
mated as 


H' om = U(ty H om U(t) 2hgaTaRe{fJo(t)} 

+hg |a(f) a 1 + a(t)* a] (b^ + bo) . (B9) 

The first term is a time dependent shift of the optical fre¬ 
quency, while the second one accounts for the linearized op¬ 
tomechanical interaction. We note that the zero-frequency 
component of the coefficients /3o(t), that we indicate with the 
symbol /3 ® 1c , and that contributes to a constant shift of the cav¬ 
ity resonance, can be taken into account by defining a renor¬ 
malized cavity frequency 

aj'=co c + 2gRe($ c } . (BIO) 

In particular we set the central frequency of the driving field 
a>i to be resonant with this shifted cavity frequency 

cj l = uj' c . fB 11) 

Correspondingly we define the time-dependent part of /3o(t) as 

o C ■ (B12) 

Explicit solutions for the steady state of a(t ) and Pj(t) can 
not be found in general, however, when a>o is sufficiently 
large JUl 


k, g\& ± \lcJa <s w 0 (B13) 

it is a valid approximation to expand the solution for a(t) and 
/ 3j(t ) at the lowest relevant order in g, and simple analytical 
expressions can be found. In detail 

a(t) = a + e 1Wof + a^e~ 1Wof + o(g 2 ) (B14) 


with 


a± 


-i£± 

k + i wq 


(B15) 


and 


A)(0 - ^o e2iWof + P^~ 2i0,0 ‘ + o(g 3 ) (B16) 


with 


Po = a ± a % v o + i (1 ± 2) cuo 1 ] _1 v 0 . (B17) 

Finally, in the interaction picture with respect to the Hamil¬ 
tonian Hq = Tiu>o Y, N j=-N^)bj, the transformed density matrix 
p — Q lH o , l^p’Q-M 0 ti'h f u ]jfjj s a mas ter equation of the form 

p = -'-\H + Hit), p\ + (£ a + £ diss ) p, (B18) 

where the time independent part of the system Hamiltonian 
H is equal to Eq. (f2Tb with the linearized coupling strength 
defined as 


G± - g a ± , 


(B19) 


while the residual time-dependent part is given by 
H(t) = 2f!gRe{y8o(f)J £a 

f=±i j =1 

+ha t [G-e- 2io>o, b 0 + G+e 21 " 0 '^ 1 ,] 

+H a [G*_e 2iwo, bl + G* + eT 2io>ot bo\ ■ (B20) 

The time-independent part is the one responsible for the 
squeezing dynamics as discussed in the main text, while H(t) 
accounts for spurious processes which tend to degrade the 
squeezing dynamics. When, consistently with Eq. (1B13I) . 
ig ± i,s|&(o|,M <s wo, the time dependent term can be ne¬ 
glected and we find Eq. ( [20l> . 


2. Optomechanics II: Simulation with an array of optical 
cavities 

A second possible implementation, complementary to the 
one just discussed, consists in the realization of the chain in 
terms of optical cavities, while the squeezed reservoir is engi¬ 
neered by optomechanical interactions between an additional 
mechanical resonator and the central cavity which in turn is 
driven by a bichromatic field as in [28]. Differently from the 
conditions identified above, here the steady state entangled 
dynamics is observed in a different optomechanical parameter 
regime corresponding to large mechanical dissipation. 

Analogously to the derivation of the previous section, we 
find that also this system can be described by a model similar 
to the one introduced in Eq. (l20ll 

P = ~\H,p] + \£^+-C ( ^: 0) ]p, (B21) 

with H formally equal to Eq. (OH . but where now bj and b\ 
are the operators for the cavity fields realizing the chain of 
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oscillators, while a and a'' are the operators for the mechanical 
resonator. Furthermore the linearized coupling strength G ± 
have a different dependence on the system parameters and are 
defined as 

G_ = gp 0 , G+ = gft* , (B22) 

with 

Po = -i £>± Vo VW ± i a> m l] -1 v 0 (B23) 

where Vo and 'W are defined in Eqs. (lB6b and (IB7 b respec¬ 
tively. We note that also the dissipative parts are slightly mod¬ 
ified. The term is equal to Eq. (0 with njj = 0, 

and describes the disspation of the optical cavities at rates jj, 
while = K(n a + 1 )D[a,cP] + Kn n D[a\a\ describes the 


dissipation of the mechanical resonator at rate k in a thermal 
bath with n a average excitations. 

In this case if the mechanical dissipation k is sufficiently 
large then the mechanical resonator can be adiabatically elim¬ 
inated and the resulting dynamics is described, by a model 
equal to Eq. ©, but where now 

r = \G+\ 2 - |G_| 2 

K 

|G_| 2 

n = n a + (2 n a + V) -y-y 

\G+\ 2 - |G_| 2 

rn -- (2n - + >) WT^J- ,B24) 
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